Synchronization of coupled demographic oscillators 
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Demographic oscillators are individual-based systems exhibiting temporal cycles sustained by the 
stochastic dynamics of the microscopic interacting particles. We here use the example of coupled 
predator-prey oscillators to show that synchronization to a common frequency can occur between 
two such systems, even if they oscillate at different frequencies in the absence of coupling. The power 
spectra of the separate and the coupled systems are computed within a van Kampen expansion in 
the inverse system size, and it is found that they exhibit two peaks at separate frequencies at low 
coupling, but that only one peak is present at large enough coupling strength. We further make 
predictions on the time behaviour of the phases of the two oscillators, and their phase difference, 
and confirm the frequency entrainment at sufficiently large coupling. Theoretical results are verified 
convincingly in numerical simulations. 
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Coupled systems of self-sustained or driven oscillators 
are abundant in physics, chemistry, biology and other 
adjacent sciences [l], and their collective behavior is of 
great interest. One of the most studied phenomena is 
that of synchronization, i.e. the adjustment of rhythms 
of oscillating units in the presence of weak interaction. 
Even if each oscillator in separation may have their own 
eigenfrequency, the ensemble of units starts oscillating at 
common frequency provided the coupling between them 
exceeds a threshold value. This phenomenon known as 
frequency entrainment or phase locking pj] . 

Synchronization effects of this type have been observed 
in a variety of contexts, going back to Christiaan Huy- 
gens in 1673 P, 0, and ranging from coupled pendula 
and clocks, to Josephson junctions, neural activity, pace- 
maker cells in the sino-atrial node to predator-prey cycles 
(see [3 and references therein). 

A variety of different mathematical models of oscil- 
lating units, and of the interactions between them has 
been considered in the literature, including coupled non- 
linear oscillators described by deterministic differential 
equations as e.g. in the celebrated Kuramoto model [1], 
time-delayed differential equations [5|], discrete systems 
0, chaotic oscillators 0], stochastic oscillatory systems 
driven by or subject to external noise [tI, Isj, and more re- 
cently individual-based spatia,l models in which diffusion 
plays the role of interaction |9[. 

In the present work we will focus on synchronization 
of oscillators driven by so-called demographic stochastic- 
ity [l3|. These are individual-based models, in which the 
constituting reagents (the individuals) interact by a set of 
simple stochastic processes resulting in the creation and 
removal of individuals, or in conversion of an individual of 
one type into an individual of another type. Chemical re- 
actions are an example of such systems, but the dynamics 



of agents in the context of game theory or the interac- 
tions in predator-prey systems T7| can be considered in 
the same framework. Depending on model parameters 
the mean-field dynamics of such systems may or may not 
allow for periodic solutions. Assuming a regime in which 
the deterministic infinite system does not exhibit oscilla- 
tions, cyclic behavior can still be found in a population 
of a finite number of individuals operating at the same 
model parameters [iTj. The random nature of the in- 
teraction between individuals is a source of stochasticity 
(so-called demographic stochasticity or intrinsic noise), 
and may under suitable circumstances lead to coherent 
sustained oscillatory behavior via a mechanism of reso- 
nant amplification. These effects have been observed in 
a number of different individual based models [l^ . and 
we will refer to such systems as 'demographic oscillators' 
in the following. 

Specifically, we will consider two coupled predator-prey 
systems, as studied in isolation in [Tl|. The two sub- 
systems will be labelled by i = 1,2, each containing two 
types of individuals, predators and prey. We will label 
predators in system z = 1 by Ai and predators in system 
i — 2 by A2, and similarly we will use Bi and B2 for the 
prey. In addition both systems may contain vacancies 
El and E2 ■ We will assume that each system contains N 
individuals (including vacancies), i.e. that the number 
of Ai, Bi and vacancies Ei sum to TV, and similarly for 
the second system. Following jTl!| the following processes 
define the dynamics of each of the sub-systems: 
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(1) 



The {(3i,ai,ji,i^i, m} are here rate constants {i = 1,2). 
The first reaction describes a birth process, the second 



and third reactions are death processes, and the remain- 
ing two reactions in each system represent predator-prey 
interaction. Up to this point the two systems i = 1,2 are 
entirely separate and it was shown in T7| that they ex- 
hibit demographic oscillations provided the reaction rates 
are chosen appropriately. In particular the power spec- 
tra of these oscillations was obtained analytically, and 
as shown in they exhibit a peak at a characteristic 
frequency set by the reaction rates. 

We will in the following focus on the case in which the 
reaction rates of the two oscillators are chosen such that 
in separation both systems oscillate at distinct frequen- 
cies. A coupling between the two systems is the intro- 
duced by allowing for the exchange of prey individuals 
according to 



Bi + E2 — > B2 + El, 
B2 + El Bi + E2- 



(2) 



e hence defines the strength of the coupling between the 
two oscillators. It is worth pointing out that the num- 
ber of individuals (including vacancies) in each system is 
constant in time under the reaction dynamics. 

For later convenience, we will denote the state of the 
coupled system by the vector n = (ni,n2, 713,714) in the 
following. The tt-q, a = 1, . . . ,4 are integers, and ni de- 
scribes the number of predators in the first system, n2 
the number of prey in the first system, and 713 and re- 
fer to predators and prey in system two. The dynamics is 
inherently stochastic and the time-evolution is governed 
by a master equation of the form 



dP(n, t) 
dt 



^[r(n|n')P(n',t)-T(n'|n)P(n,t)], (3) 



where r(n|n') is the transition rate from state n' to n. 
These rates are defined by the set of allowed reactions, 
and we will not report their detailed mathematical form. 

Similar to the mean field limit results in de- 

terministic equations for the mean values xi = ni/N, 
Ui = n2/N , X2 = ns/N and 1/2 — n^/N. For the coupled 
system these equations are of the form 

ii = 2z^ia;i?/i — aiXi 

yi = 2l3i{l - xi - yi)yi - ^lyi - 2{vi + iii)xiyi 
+e [2/2(1 -xi- yi) - yi(l - X2 - 2/2)] 

±2 = 2/iiX2j/2 - 722^2 

2)2 = 2/32(1 - 2^2 - 2/2)2/2 - 72?/2 - 2(1/2 + ^^2)2:22/2 

+e [2/1(1 - a;2 - 2/2) - 2/2(1 - a;i - 2/1)] • (4) 

As the only coupling between the two predator-prey sys- 
tem is through the exchange of prey individuals, the 
equations for the pairs {xi,yi) and (x2, 2/2) respectively 
decouple for e = 0. As discussed in [11[ deterministic 
systems of this type approach a non-trivial fixed point 
asymptotically at £ = and do not allow any limit-cycle 
solutions. 
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FIG. 1: Predator densities in the two subsystems as a func- 
tion of time. Upper two panels: the two oscillators are cou- 
pled and synchronize; lower two panels: uncoupled case, no 
synchronization. Horizontal lines in each panel represent the 
fixed-point solution of the deterministic dynamics Eqs. Q, 
noisy curves show a single simulation run of the stochastic 
system with N = 5000 individuals in each oscillator. 



The above mean-field equations are however valid only 
in the limit of infinite systems, N ~^ 00. At finite sizes 
coherent amplified oscillations about the mean-field fixed 
points have been observed pdj . and their power spectra 
have been computed within a systematic expansion in 
the inverse system size. The spectra exhibit peaks at a 
non-zero frequency, determined by the rate constants and 
setting the characteristic period of the stochastic oscilla- 
tions. We will focus on cases in which the two oscilla- 
tors i = 1,2 have characteristically different frequencies 
of their respective demographic oscillations. Specifically 
we will use (ai, 71, i/i, ^ii) = (0.1,0.1,0.0,0.25,0.05) 
and (a2,/32,72, i/2,/i2) = (0.1,0.1,0.1,0.5,0.5). In sepa- 
ration the two systems show cyclic behavior as shown 
in the lower two panels of Fig. [TJ The periodic be- 
havior of the system at finite sizes can be investigated 
analytically within a system-size expansion as proposed 
first by van Kampen In essence the method con- 

sists in expanding about the mean-field solution, writing 
e.g. ni/N = xi + zi/^/N and similarly for the remain- 
ing three components of the system. Changing from the 
variables 711, 712, ^3, to zi, Z2, Z3, Z4 in the master equa- 
tion one then performs a systematic expansion in powers 
of N^^/"^. We will not present the details of the mathe- 
matics here, as these have been reported for similar sys- 
tems in the literature [ll|, [l3j llHl ; but will only report 
the final result. The leading order of this expansion re- 
produces the mean-field equations as expected. In 
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FIG. 2: (color on-line) Power spectra Pii(lj) (blue circles) 
and P33 {(jj) (red squares) of the predator densities in the two 
sub-systems at diflterent magnitudes of the coupling strength 
e. Solid lines are obtained from the theory, markers from 
simulations of the individual-based model (A^ = 10'*, averages 
over at least 200 samples are taken in each panel). 



next-to-leading order one obtains a Langevin equation 



z = Jz + r]{t) 



(5) 



for the vector z — (21,22,-^3,^4) describing fluctuations 
about the mean-field result, rj is here a 4-component 
Gaussian noise vector, which comes out as white in time, 
but with correlations between its components, as we will 
specify below. Focusing on the asymptotic regime in 
which the mean-field theory attains a fixed point, J is 
the 4x4 Jacobian of the system ^ at this fixed point 
[m . The correlator of the noise variables can be worked 
out and one finds 



12 



{Vait)r]b{t')) = S{t - t') 2J ariSna)r{Snb) 



(6) 



for the case of two coupled predator-prey systems, a, 6 G 
{1, . . . , 4} label the components of the noise, and the in- 
dex r runs through the twelve possible reactions in the 
system (five in each subsystem, and the two exchange 
reactions), stands for the rate at which reaction r 
occurs in the mean-field system. The integer numbers 
{dna)r finally denote the number of particles of type 
a that are created {{Sna)r > 0), or removed from the 
system {{dna)r < 0) during one occurrence of reaction 

r e {1, . . . , 12}. If we give reaction Bi + Ei ^ Bi + Bi 
the label r = 1 then one has ai ~ /3ij/i(l — xi — j/i), 
and {dni)i — 0, {Sn2)i — 1, (^^13)1 = (<5n4)i = for ex- 
ample. Analogous expressions can be formulated for the 
remaining reactions. The above Langevin equation for z 
is linear, and can be solved in Fourier space. Details for 



similar systems are again found in the literature, so that 
we do not describe the details of the intermediate steps 
here. The outcome is the set of power spectra 

Pab{uj) ^ {Za{~Uj)%{uj)) , fl, 6 G { 1 , . . . , 4} , (7) 

where Za{uj) is the Fourier transform of Za{t), and 
where (• • •) stands for an average over realizations of the 
Langevin dynamics. Results from this analytical com- 
putation are shown in Fig. [2l where we plot the power 
spectra of the predator densities in either of the two sub- 
systems, Pii(w) and P33(a;). The solid fines are the pre- 
dictions of the theory, markers are from simulations based 
on the Gillespie algorithm [l^l- As seen in the figure the 
agreement between theory and numerical experiment is 
excellent. With our choice of parameters the dominating 
frequency components in the fluctuations of the number 
of predators are quite distinct in the absence of coupling 
between the systems. Two peaks at distinct frequencies 
are seen in the top left panel of Fig. O As we gradually 
increase the coupling e between the two demographic os- 
cillators, the positions of the two peaks approach each 
other, and flnally the dominant frequencies of both sys- 
tems essentially coincide at large enough e (see Fig. 
leading to frequency entrainment. This effect is also vis- 
ible in the time series shown in Fig. [T] While the fluctu- 
ations in the number of predators have characteristically 
distinct periods in the absence of coupling between the 
two oscillators (lower panel), coherent oscillations at a 
unique common frequency are found for e = 0.1 (upper 
panel). Due to the stochasticity in the oscillations this 
coherence can never be expected to be perfect, but the 
qualitative effect is certainly visible in the figure. 

To further illustrate the synchronization effect in the 
system we have monitored the position of the peaks in 
the power spectra (as obtained from the system-size ex- 
pansion) as s is varied. Results are shown in the upper 
panel of Fig. [3l and while the peaks are well separated in 
the regime of small coupling, their peak positions essen- 
tially coincide for values of the coupling strength greater 
than about e = 0.07. While this is not an analytically 
exact result, the deviation between the position of the 
peaks is minute in this regime, and for all practical pur- 
poses it is therefore fair to conclude that the two systems 
have indeed synchronized. 

Further insight can be gained by computing the eigen- 
values of the Jacobian J at the fixed-point of the de- 
terministic system. The presence of sustained demo- 
graphic oscillations in the stochastic system is here in- 
dicated by non-zero imaginary parts of these eigenval- 
ues. In essence, non- vanishing imaginary parts cause an 
oscillatory approach of the deterministic system to the 
stable fixed point (real parts of the eigenvalues are neg- 
ative). In the finite system the stochastic nature of the 
microscopic dynamics results in an ongoing random per- 
turbation away from this fixed point. An instantaneous 
perturbation would decay in an oscillatory manner back 
to the fixed point, but as perturbations due to the de- 
mographic noise occur at all times, a coherent oscillatory 



4 



0.14 
0.12 

o 
CI. 

3 0.1 
0.08 



1 ' 1 
system 1 


- 


system 2 




1 , 1 





"0 



0.05 



0.1 






0.05 


0.1 


1 , 1 , 














I.I. 



FIG. 3: (color on-line) Position of the peaks of the power 
spectra of predator fluctuations (upper panel) as a function of 
the coupling strength e. The lower panel shows the imaginary 
parts of the eigenvalues of the Jacobian J (see text for details). 




FIG. 4: (color on-line) Instantaneous phases of the preda- 
tor fluctuations in the two oscillators. Symbols are from 
simulations (circles are (^i, squares cj>3, diamonds — cj>3\ 
N — 10^, averaged over 20 samples), solid lines have slopes 

= CLlo'lo' Paa{l^')/[j^ dw'Paa (w')] , wherC Paa{u) (o £ 

{1,3}) is the analytically computed spectrum, n « 1.2 is 
a cutoff used to reduce the effects of phase slips [l^. Mea- 
surements of the phase start ai t = 500, to allowing for equi- 
libration. 



pattern emerges. In a crude approximation one may es- 
timate that the frequency of such oscillations is set by 
the imaginary parts of the eigenvalues of the Jacobian. 
Our system is 4-dimensional, so we have four eigenval- 
ues of the Jacobian of the coupled system, and three 
different scenarios are possible: (i) all four eigenvalues 
are real (no demographic oscillations), (ii) they form two 
distinct pairs of complex conjugates (resulting in two dis- 
tinct characteristic frequencies), or (iii) two of them are 
real, and the remaining two are complex conjugates with 
non-zero imaginary part (resulting in only one charac- 
teristic frequency). As seen in the lower panel of Fig. 
131 (ii) is the case at low coupling strengths, up to about 
£ = 0.11 for the specific reaction rates chosen throughout 
this paper. Above this critical coupling strength, only 
one pair of complex-conjugate eigenvalues is found, indi- 
cating that there is only one characteristic frequency in 
the system. In case of two pairs of complex-conjugates, 
care needs to be taken however in identifying the result- 
ing frequencies with those of oscillations of the various 
species concentrations in the system. The precise eigen- 
vector structure may here be important if one wanted to 
identify the relevant normal modes (no such attempt has 
been made here), and we must stress that the peak of 
the different power spectra is generally not located pre- 
cisely at the frequencies set by the imaginary parts of the 
eigenvalues, but that their real parts and the correlations 
between the noise components 77^ may have an impact on 
their position as well. 

To further study the phase synchronization of the two 
demographic oscillators we have measured their relative 
phases numerically, and report results in Fig. [H The 



phase (/)(t) of a narrow frequency-band periodic signal is 
here defined so that it increases by 27r within one oscil- 
lation cycle, and grows in proportion to the fraction of 
the period that has elapsed |lj|. Extracting the instan- 
taneous phase 0(t) from a stochastic signal generated by 
Gillespie simulations is non-trivial, we have here resorted 
to Hilbert transform techniques as described in [l| , note 
also [IBl . As seen in the left panel of Fig. [Hthe analysis of 
the phases of the two demographic systems confirms that 
they oscillate at distinctively different frequencies if un- 
coupled, and the difference in their phases increases lin- 
early in time. At sufficiently large coupling (right panel) 
the phases of both sub-systems grow at the same slope 
in time (the slopes of the curves in the right panel of Fig. 
2]) coincide within a margin of about 2 per cent) , and the 
phase difference hence remains finite in time, indicating 
phase locking and frequency entrainment. 

Before concluding, we would like to point out that a 
systematic investigation of the behaviour of the coupled 
predator-prey system at values of the reaction rates dif- 
ferent from the choice made in the present work is still 
pending, and no claim is made here that synchronization 
will occur for all choices of the reaction rates in In- 
deed as shown in Fig. [5] oscillation death 1(5] may occur 
at sufficiently large coupling strength, i.e. the (global) 
maximum of the power spectra of individual species may 
no longer be located at a non-zero frequency, indicating 
that demographic oscillations no longer exists for certain 
species in the system. For other choices of the reaction 
rates, this may in principle occur before synchronization 
is reached, so that such a system would not exhibit phase 
locking, but rather oscillation death would preempt syn- 
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FIG. 5: (color on-line) Power spectra of the fluctuations of 
the four species at £ = 0.3. Solid lines are analytical results 
obtained within the van Kampen expansion, markers from 
simulations at TV = 10* (averaged over 500 samples). 



chronization of demographic oscillations under such cir- 



cumstances. 

In summary we have shown that the paradigm of syn- 
chronization at sufficiently large coupling extends to de- 
mographic oscillators, i.e. system in which the periodic 
behavior in each unit is sustained by coherent amplifica- 
tion of demographic noise. The model system we have 
focused on is that of two coupled predator-prey systems, 
and the synchronization effects we find indicate that the 
entrained dynamics of lynx abundances in different re- 
gions of Canada as reported in [l^ may be explained 
within the picture of quasi-cycles driven by intrinsic noise 
(see also [19] and references therein for theoretical models 
of synchronization in spatially extended ecological sys- 
tems) . We expect that a similar mechanism will occur in 
other demographic oscillators, for example in biochemical 
reactions [l3] or in models of epidemics [l2| , and in sys- 
tems which are composed of more than two self-sustained 
oscillators. 
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